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Abstract 

This  article  provides  an  overview  of  the  effort  centered  on  the  HYbrid  Coordinate  Ocean  Model  (HYCOM)  to  develop  an  eddy- 
resolving,  real-time  global  and  basin-scale  ocean  hindcast,  nowcast,  and  prediction  system  in  the  context  of  the  Global  Ocean  Data 
Assimilation  Experiment  (GODAE),  The  main  characteristics  of  HYCOM  are  first  presented,  followed  by  a  description  and 
assessment  of  the  present  near  real-time  Atlantic  forecasting  system.  Regional/coastal  applications  are  also  discussed  since  an 
important  attribute  of  the  data  assimilative  HYCOM  simulations  is  the  capability  to  provide  boundary  conditions  to  regional  and 
coastal  models.  The  final  section  describes  the  steps  taken  toward  the  establishment  of  the  fully  global  eddy-resolving  HYCOM 
data  assimilative  system  and  discusses  some  of  the  difficulties  associated  with  advanced  data  assimilation  given  the  size  of  the 
problem. 

©  2006  Elsevier  B,V«  All  rights  reserved. 
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1.  Introduction 

A  broad  partnership  of  institutions1  is  presently 
collaborating  in  developing  and  demonstrating  the 
performance  and  application  of  eddy-resolving,  real¬ 
time  global  and  basin-scale  ocean  hindcast,  nowcast,  and 
prediction  systems  using  the  HYbrid  Coordinate  Ocean 

*  Corresponding  author. 

E-mail  address:  cchasstgnet@coaps.fsu.edu  (El  Chassignet), 

1  U,  of  Miami*  NRL,  Los  Alamos,  NQAA/NCEP,  NOAA7AOML, 
NOAA/PMEL,  PSJ,  FNMOC,  NAVOCEANO,  SHOM,  LEG1,  OPeN- 
DAP,  U.  of  North  Carolina,  Rutgers,  U.  of  South  Florida,  Fugro- 
GEOS,  ROFFS,  Orbimage,  Shell,  ExxonMobil 

0924 -7 963/S  -  see  front  matter  ©  2006  Elsevier  B.V,  All  rights  reserved, 
doi ;  1 0. 1 01 6/j  j  marsy  s*2O05 .09. 0 1 6 


Model  (HYCOM),  The  plan  is  to  transition  these  systems 
for  operational  use  by  the  U,S,  Navy  at  the  Naval 
Oceanographic  Office  (NAVOCEANO),  Stennis  Space 
Center,  MS,  and  the  Fleet  Numerical  Meteorology  and 
Oceanography  Center  (FNMOC),  Monterey,  CA;  and 
by  NOAA  at  the  National  Centers  for  Environmental 
Prediction  (NCEP),  Washington,  D.C  The  partnership  is 
also  the  eddy-resolving  global  ocean  data  assimilative 
system  development  effort  that  is  sponsored  by  the  US. 
component  of  the  Global  Ocean  Data  Assimilation 
Experiment  (GODAE).  GODAE  is  a  coordinated  inter¬ 
national  effort  envisioning  4ia  global  system  of  observa¬ 
tions ,  communications .  modeling,  and  assimilation  that 
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will  deliver  regular  comprehensive  information  on  the 
state  of  the  oceans,  in  a  way  that  will  promote  and 
engender  wide  utility  and  availability  of  this  resource  for 
maximum  benefit  to  the  community  9\  Three  of  the 
GODAE  specific  objectives  are  to  apply  state-of-the-art 
models  and  assimilation  methods  to  produce  short-range 
open  ocean  forecasts,  boundary  conditions  to  extend 
predictability  of  coastal  and  regional  subsystems,  and 
initial  conditions  for  climate  forecast  models  (GODAE 
Strategic  Plan,  International  GODAE  Steering  Team, 
2000).  HYCOM  development  is  the  result  of  collabora¬ 
tive  efforts  among  the  University  of  Miami,  the  Naval 
Research  Laboratory  (NRL),  and  the  Los  Alamos 
National  Laboratory  (LANL),  as  part  of  the  multi- 
institutional  HYCOM  Consortium  for  Data  Assimilative 
Ocean  Modeling  funded  by  the  National  Ocean 
Partnership  Program  (NOPP)  in  1999  to  develop  and 
evaluate  a  data  assimilative  hybrid  isopycnabsigma- 
pressure  (generalized)  coordinate  ocean  model  (Bleck, 
2002;  Chassignet  et  ah,  2003  ;  Hal  li well,  2004), 

Numerical  modeling  studies  over  the  past  several 
decades  have  demonstrated  advances  in  both  model 
architecture  and  the  availability  of  computational  re¬ 
sources  for  the  scientific  community.  Perhaps  the  most 
noticeable  aspect  of  this  progression  has  been  the 
evolution  from  simulations  on  coarse-resolution  hori¬ 
zontal/vertical  grids  outlining  basins  of  simplified  geo¬ 
metry  and  bathymetry  and  forced  by  idealized  stresses, 
to  fine-resolution  simulations  incorporating  realistic 
coastal  definition  and  bottom  topography,  forced 
by  observational  data  on  relatively  short  time  scales 
(Hurlburt  and  Hogan,  2000;  Smith  et  al.,  2000; 
Chassignet  and  Garraffo,  2001).  Traditional  Ocean 
General  Circulation  Models  (OGCMs)  use  a  single  co¬ 
ordinate  type  to  represent  the  vertical,  but  recent  model 
comparison  exercises  performed  in  Europe  (DYnamics 
of  North  Atlantic  MOdels— DYNAMO)  (Willebrand 
et  al.,  2001)  and  in  the  U.S.  (Data  Assimilation  and 
Model  Evaluation  Experiment— DAM^E)  (Chassignet 
et  al.,  2000)  have  shown  that  no  single  vertical  coor¬ 
dinate-depth,  density,  or  terrain-following  cr-levels-can 
by  itself  be  optimal  everywhere  in  the  ocean.  These  and 
earlier  comparison  studies  (Chassignet  et  ah,  1996; 
Roberts  et  al.,  1996;  Marsh  et  al.,  1996)  have  shown  that 
the  models  considered  are  able  to  simulate  the  large-scale 
characteristics  of  the  oceanic  circulation  reasonably  well, 
but  that  the  interior  water  mass  distribution  and  asso¬ 
ciated  thermohaline  circulation  are  strongly  influenced 
by  localized  processes  that  are  not  represented  equally  by 
each  model's  vertical  discretization.  The  choice  of  the 
vertical  coordinate  system  is  one  of  the  most  important 
aspects  of  an  ocean  model's  design  and  practical  issues 


of  representation  and  parameterization  are  often  directly 
linked  to  the  vertical  coordinate  choice  (Griffies  et  al., 
2000),  Currently,  there  are  three  main  vertical  coordi¬ 
nates  in  use,  none  of  which  provides  universal  utility. 
Hence,  many  developers  have  been  motivated  to  pursue 
research  into  hybrid  approaches.  Isopycnal  (density 
tracking)  layers  are  best  in  the  deep  stratified  ocean, 
z-Ievels  (constant  fixed  depths)  are  best  used  to 
provide  high  vertical  resolution  near  the  surface  with¬ 
in  the  mixed  layer,  and  ex-levels  are  often  the  best 
choice  in  shallow  coastal  regions.  HYCOM  combines 
all  three  approaches  and  the  optimal  distribution  is 
chosen  at  every  time  step.  The  model  makes  a  dy¬ 
namically  smooth  transition  between  the  coordinate 
types  via  the  continuity  equation  using  the  hybrid 
vertical  coordinate  generator. 

The  layout  of  the  paper  is  as  follows.  First,  in  Section  2, 
an  overview  of  the  HYCOM  characteristics  is  presented 
with  model  performance  illustrated  using  non-data 
assimilative  basin-scale  and  regional  nested  simulations. 
The  near  real-time  North  Atlantic  Ocean  data  assimilative 
system  in  then  introduced  in  Section  3  and  its  hindcast 
capabilities  evaluated.  In  Section  4,  issues  associated  with 
regional/coastal  applications  are  introduced  and  discussed. 
Future  development  plans  are  presented  in  Section  5. 

2.  The  ocean  model 

HYCOM  is  designed  to  provide  a  significant  im¬ 
provement  over  existing  operational  OGCMs,  since  it 
overcomes  design  limitations  of  present  systems  as  well 
as  limitations  in  vertical  discretization.  The  ultimate  goal 
is  a  more  streamlined  system  with  improved  perfor¬ 
mance  and  an  extended  range  of  applicability  (e.g.,  the 
present  U.S.  NAVY  systems  are  seriously  limited  in 
shallow  water  and  in  handling  the  transition  from  deep  to 
shallow  water).  The  generalized  coordinate  (hybrid) 
ocean  model  HYCOM  used  in  this  study  retains  many 
of  the  characteristics  of  its  predecessor,  the  isopycnic 
coordinate  model  MICOM  (Miami  Isopycnic  Coordi¬ 
nate  Model)  (Bleck  et  al.,  1992;  Bleck  and  Chassignet, 
1994),  while  allowing  coordinate  surfaces  to  locally 
deviate  from  isopycnals  wherever  the  latter  may  fold, 
outcrop,  or  generally  provide  inadequate  vertical  reso¬ 
lution  in  portions  of  the  model  domain.  Hybrid  coor¬ 
dinates  can  mean  different  things  to  different  people:  it 
can  be  a  linear  combination  of  two  or  more  conventional 
coordinates  (Song  and  Haidvogel,  1994;  Ezer  and 
Mellor,  2004;  Barron  et  al.,  in  press)  or  truly  generalized, 
i.e.  aims  to  mimic  different  types  of  coordinates  in 
different  parts  of  a  model  (Bleck,  2002;  Burchard  and 
Beckers,  2004;  Adcroft  and  Hallberg,  2006;  Song  and 
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Hou,  2006),  HYCGM  uses  the  same  equations  as 
MICOM>  except  that  they  have  been  modified  to  account 
for  nonzero  horizontal  density  gradient  within  all  layers, 
not  just  the  top  layer  as  in  MICGM.  HYCOM  remains  a 
Lagrangian  layer  model  in  the  sense  that  the  MICOM 
solution  procedure  is  unmodified,  except  that  remapping 
of  the  vertical  coordinate  is  performed  via  a  hybrid 
coordinate  generator  at  the  end  of  each  baroclinic  time 
step  (Bleck,  2002;  Halliwell,  2004).  HYCOM  is  thus 
classified  as  a  Lagrangian  Vertical  Direction  (LVD) 
model  where  the  continuity  (thickness  tendency)  equa¬ 
tion  is  solved  prognostically  throughout  the  domain, 
with  the  Arbitrary  Lagrangian-Eulerian  (ALE)  technique 
used  to  remap  the  vertical  coordinate  and  maintain  dif¬ 
ferent  coordinate  types  within  the  domain  (Adcroft  and 
Hallberg,  2006).  This  differs  from  Eulerian  Vertical 
Direction  (E  VD)  models  with  fixed  z-  and  cr-coordinates 
that  use  the  continuity  equation  to  diagnose  vertical 
velocity. 

The  freedom  to  adjust  the  vertical  spacing  of  the 
coordinate  surfaces  in  HYCOM  simplifies  the  numerical 
implementation  of  several  physical  processes  (mixed  layer 
detrainment,  convective  adjustment,  sea  ice  modeling,  „.) 
without  robbing  the  model  of  the  basic  and  numerically 
efficient  resolution  of  the  vertical  that  is  characteristic  of 
isopycnic  models  throughout  most  of  the  ocean's  volume 
(see  Section  2 A  for  details).  The  capability  of  assigning 
additional  coordinate  surfaces  to  the  oceanic  mixed  layer 
in  HYCOM  allows  the  option  of  implementing  sophisti¬ 
cated  vertical  mixing  turbulence  closure  schemes.  The 
latest  release  of  HYCOM  has  five  primary  vertical  mixing 
algorithms,  of  which  three  are  vertical  diffusion  models 
and  two  are  slab  models  (see  Section  22  for  details).  The 
choice  of  the  vertical  mixing  parameterization  is  also  of 
importance  in  areas  of  strong  entrainment,  such  as 
overflows  (see  Section  2,3  for  details). 

2.  L  Hybrid  coordinate  generator  and  its  transition  to 
coastal  regions 

The  implementation  of  the  generalized  vertical  co¬ 
ordinate  in  HYCOM  follows  the  theoretical  foundation 
set  forth  in  Bleck  and  Boudra  (1981)  and  Bleck  and 
Benjamin  (1993):  i.e.,  each  coordinate  surface  is 
assigned  a  reference  isopycnal.  The  model  continually 
checks  whether  or  not  grid  points  lie  on  their  reference 
isopycnals  and,  if  not,  attempts  to  move  them  vertically 
toward  the  reference  position.  However,  the  grid  points 
are  not  allowed  to  migrate  when  this  would  lead  to 
excessive  crowding  of  coordinate  surfaces.  Thus, 
vertical  grid  points  can  be  geometrically  constrained  to 
remain  at  a  fixed  depth  while  being  allowed  to  join  and 


follow  their  reference  isopycnals  in  adjacent  areas 
(Bleck,  2002).  The  default  configuration  in  HYCOM 
is  one  that  is  isopycnal  in  the  open  stratified  ocean,  but 
makes  a  dynamically  smooth  transition  to  a  coordinates 
in  shallow  coastal  regions  and  to  fixed  pressure-level 
coordinates  (hereafter  referred  to  as  p)  in  the  surface 
mixed  layer  and/or  unstratified  seas  (Fig,  1 ).  In  doing  so, 
the  model  combines  the  advantages  of  the  different 
coordinate  types  in  optimally  simulating  coastal  and 
open-ocean  circulation  features.  It  is  left  to  the  user  to 
define  the  coordinate  separation  constraints  that  control 
regional  transitions  among  the  three  coordinate  choices 
as  described  in  the  appendix. 

After  the  model  equations  are  solved,  the  hybrid 
coordinate  generator  then  relocates  vertical  interfaces  to 
restore  isopycnic  conditions  in  the  ocean  interior  to  the 
greatest  extent  possible  while  enforcing  the  minimum 
thickness  requirements  specified  by  ( 1 )  in  the  appendix.  If 
a  layer  is  less  dense  than  its  isopycnic  reference  density, 
the  generator  attempts  to  move  the  bottom  interface 
downward  so  that  the  flux  of  denser  water  across  this 
interface  increases  density.  If  the  layer  is  denser  than  its 
isopycnic  reference  density,  the  generator  attempts  to 
move  the  upper  interface  upward  to  decrease  density.  In 
both  cases,  the  generator  first  calculates  the  vertical 
distance  over  which  the  interface  must  be  relocated  so  that 
volume- weighted  density  of  the  original  plus  new  water  in 
the  layer  equals  the  reference  density.  The  minimum 
permitted  thickness  of  each  layer  at  each  model  grid  point 
is  then  calculated  using  (1)  in  the  appendix.  The  final 
minimum  thickness  is  then  calculated  using  a  “cushion” 
function  (Bleck,  2002)  that  produces  a  smooth  transition 
from  the  isopycnic  to  the p  and  a  domains.  The  minimum 
thickness  constraint  is  not  enforced  at  the  bottom  in  the 
open  ocean,  permitting  the  model  layers  to  collapse  to 
zero  thickness  there,  as  in  MICOM.  Repeated  execution 
of  this  algorithm  at  every  time  step  maintains  layer  density 
very  close  to  its  reference  value  as  long  as  a  minimum 
thickness  does  not  have  to  be  maintained  and  diabatic 
processes  are  weak.  To  insure  that  a  permanent  p - 
coordinate  domain  exists  near  the  surface  year  round  at 
all  model  grid  points,  the  reference  densities  of  the 
uppermost  layers  are  assigned  values  smaller  than  any 
density  values  found  in  the  model  domain. 

Fig.  1  illustrates  the  transition  that  occurs  between  p/cr 
and  isopycnic  ( p )  coordinates  in  the  fall  and  spring  in  the 
upper  400  m  and  over  the  shelf  in  the  East  China  and 
Yellow  Seas.  In  the  fall,  the  water  column  is  stratified  and 
can  be  largely  represented  with  isopycnals;  in  the  spring, 
the  water  column  is  homogenized  over  the  shelf  and  is 
represented  by  a  mixture  of  p  and  a  coordinates.  A 
particular  advantage  of  isopycnic  coordinates  is  illustrated 
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Fig.  1.  Upper  400  m  north-south  velocity  cross-section  along  124.5DE  in  a  l/25d  East  China  and  Yellow  Seas  HYCOM  embedded  in  a  1/6°  North 
Pacific  configuration  forced  with  climatological  monthly  winds:  (a)  in  the  fall,  the  water  column  is  stratified  over  the  shelf  and  can  be  represented  with 
isopycnals  (p);  {b)  in  the  spring,  the  water  column  is  homogenized  over  the  shelf  and  the  vertical  coordinate  becomes  a  mixture  of  pressure  (p)  levels 
and  terrain- following  (er)  levels.  The  isopycnic  layers  arc  numbered  over  the  shelf,  the  higher  the  number,  the  denser  the  layer. 


by  the  density  front  formed  by  the  Kuroshio  above  the 
peak  of  the  sharp  (tip)  topography  at  the  shelfbreak  in 
Fig.  la.  Since  the  lip  topography  is  only  a  few  grid  points 
wide,  this  topography  and  the  associated  front  is  best 
represented  in  isopycnic  coordinates.  In  other  applications 
in  the  coastal  ocean,  it  may  be  more  desirable  to  provide 
high  resolution  from  surface  to  bottom  to  adequately 
resolve  the  vertical  structure  of  water  properties  and  of  the 
bottom  boundary  layer.  Since  vertical  coordinate  choices 
for  open  ocean  HYCOM  runs  typically  maximize  the 
fraction  of  the  water  column  that  is  isopycnic,  it  is  often 
necessary  to  add  more  layers  in  the  vertical  to  coastal 
HYCOM  simulations  nested  within  larger-scale  HYCOM 
runs.  The  nested  West  Florida  Shelf  simulations  analyzed 
by  Halliwell  (in  preparation)  use  this  technique,  which  is 
illustrated  in  the  cross-sections  in  Fig.  2.  The  original 
vertical  discretization  is  compared  to  two  others  with  six 
layers  added  at  the  top:  one  with  p  coordinates  and  the 
other  with  a  coordinates  over  the  shelf.  This  illustrates  the 
flexibility  with  which  vertical  coordinates  can  be  chosen 
using  the  minimum  layer  thickness  algorithm  in  the 
appendix,  Halliwell  (in  preparation)  documents  the 
advantages  of  using  high-resolution  a  coordinates 
compared  to  the  other  two  choices  shown  in  Fig.  2, 


Maintaining  hybrid  vertical  coordinates  can  be  thought 
of  as  upwind  finite  volume  advection.  The  original  grid 
generator  (Bleck,  2002)  used  the  simplest  possible 
scheme  of  this  type,  the  1st  order  donor-cell  upwind 
scheme,  A  major  advantage  of  this  scheme  is  that  moving 
a  layer  interface  does  not  affect  the  layer  profile  in  the 
down- wind  (detraining)  layer,  which  greatly  simplifies 
remapping  to  isopycnal  layers.  However,  the  scheme  is 
diffusive  when  layers  are  remapped  (there  is  no  diffusion 
when  layer  interfaces  remain  at  their  original  location). 
Isopycnal  layers  require  minimal  remapping  in  response 
to  weak  interior  diapycnal  diffusivity,  but  fixed  coordinate 
layers  often  require  significant  remapping,  especially  in 
regions  with  significant  upwelling  or  down  welling. 
Therefore,  to  minimize  diffusion  associated  with  the 
remapping,  the  grid  generator  now  uses  the  piecewise 
linear  method  with  a  monoton  ized  central -difference 
(MC)  limiter  (van  Leer,  1997)  for  layers  that  are  in  fixed 
coordinates  white  still  using  donor-cell  upwind  for  layers 
that  are  non- fixed  (and  hence  tending  to  isopycnal 
coordinates).  The  piecewise  linear  method  replaces  the 
“constant  within  each  layer”  profile  of  donor-cell  with  a 
linear  profile  that  equals  the  layer  average  at  the  center  of 
the  layer.  The  slope  must  be  limited  to  maintain 
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Fig,  2,  Cross-sections  of  layer  density  and  model  interfaces  across  the  West  Florida  Shelf  illustrating  the  capability  to  add  new  layers  at  the  top  for  a 
nested  coastal  simulation  and  the  capability  to  specify  different  coordinate  types  over  the  shelf  The  1/25°  West  Florida  Shelf  subdomain  covers  the 
Gulf  of  Mexico  east  of  87°W  and  north  of  23°N  and  is  embedded  in  a  1/2 5°  lntra-Amencas  Sea,  itself  nested  within  a  cl imaiologi tally- forced  1/12° 
Atlantic  basin  HYCOM  simulation  (HalliweJL  in  preparation). 


monotonicity;  there  are  many  possible  limiters  but  the  MC 
limiter  is  one  of  the  more  widely  used  (Leveque,  2002). 

2.2.  Mixed  layer  options 

As  noted  earlier,  the  capability  of  assigning  addi¬ 
tional  coordinate  surfaces  to  the  HYCOM  mixed  layer 
allows  the  option  of  implementing  sophisticated  vertical 
mixing  turbulence  closure  schemes  (see  Halliwell,  2004 
for  a  review).  The  full  set  of  vertical  mixing  options 
contained  in  the  latest  version  of  HYCOM  (http:// 
hycom.rsmas.miami.edu)  includes  five  primary  vertical 
mixing  submodels,  of  which  three  are  “continuous” 
vertical  diffusion  models  and  two  are  predominantly  or 
totally  bulk  models.  The  three  vertical  diffusion  models. 


which  govern  vertical  mixing  throughout  the  water 
column,  are  the  K- Profile  Parameterization  of  Large 
et  al.  (1994)  (KPP),  the  level  2.5  turbulence  closure  of 
Meilor  and  Yamada  (1982)  (MY),  and  the  Goddard 
Institute  for  Space  Studies  (GISS)  level  2  turbulence 
closure  of  Canuto  et  al.  (2001 , 2002).  The  other  two  are 
the  quasi-bulk  dynamical  instability  submodel  of  Price 
et  al.  (1986)  (PWP)  and  the  bulk  Kraus  and  Turner 
(1967)  submodel  (KT).  Since  these  latter  two  mixed 
layer  models  do  not  provide  mixing  from  surface  to 
bottom,  HYCOM  contains  two  diapycnal  mixing 
models,  one  explicit  and  one  implicit,  to  provide  this 
mixing  in  the  interior  ocean.  All  mixing  schemes  within 
HYCOM  are  kept  up  to  date.  The  MY  model  is  the 
version  implemented  in  the  Princeton  Ocean  Model 
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(POM),  specifically  POM98.  The  latest  recommended 
coefficients  are  implemented  in  KPP,  and  we  have  an 
ongoing  collaboration  with  NASA/GISS  to  implement 
the  latest  changes  in  their  model.  Future  plans  include 
implementing  and  testing  new  mixing  models,  such  as 
the  generic  length  scale  equation  turbulence  closure  of 
Umlauf  and  Burchard  (2003). 

The  following  procedure  is  used  to  implement  the 
three  vertical  diffusion  submodels  (KPP,  GISS,  and 
MY).  Velocity  components  are  interpolated  to  the  p 
grid  points  from  their  native  u  and  v  points.  The  one¬ 
dimensional  submodels  are  then  run  at  each  p  point  to 
calculate  profiles  of  viscosity  coefficients  along  with  T 
and  S  diffusion  coefficients  on  model  interfaces.  The 
one-dimensional  vertical  diffusion  equation  is  then 
solved  at  each  p  point  to  mix  I \  5,  and  tracer  variables, 
which  involves  the  formulation  and  solution  of  a  tri¬ 
diagonal  matrix  system  using  the  algorithm  provided 
with  the  KPP  submodel  (Large  et  al.,  1994,  1997).  To 
mix  momentum  components,  viscosity  profiles  stored 
on  interfaces  at  p  grid  points  are  horizontally  interpo¬ 
lated  to  interfaces  at  u  and  v  grid  points;  then  the  vertical 
diffusion  equation  is  solved  on  both  sets  of  points.  These 
three  mixing  models  all  diagnose  mixed  layer  thickness 
using  the  method  of  Kara  et  al.  (2000,  2003),  which  is 
implemented  as  follows:  the  user  first  specifies  a 
minimum  temperature  jump  for  estimating  mixed  layer 
thickness,  which  is  converted  to  an  equivalent  density 
jump  using  the  equation  of  state.  The  mixed  layer  base  is 
then  assumed  to  reside  at  the  depth  where  density  differs 
from  layer  1  density  by  the  value  of  this  jump.  Moving 
down  from  layer  1,  the  first  model  layer  where  the 
density  exceeds  layer  1  density  by  more  than  the  value 
of  this  jump  is  identified.  Given  the  central  depth  and 
density  of  this  layer  along  with  the  central  depth  and 
density  of  the  layer  above,  linear  interpolation  is  used  to 
estimate  the  thickness. 

The  original  KPP  submodel  did  not  contain  a  bottom 
boundary  layer  parameterization,  but  one  was  added  by 
Hal  I  i  well  (in  preparation)  to  perform  coastal  simulations 
over  the  West  Florida  Shelf  It  is  essentially  the  same 
parameterization  used  for  the  surface  layer,  but  turned 
^upside-down”.  The  procedure  implemented  in  HYCOM 
essentially  follows  the  procedure  implemented  in  the 
Regional  Ocean  Model  System  (ROMS)  (Durski  et  al., 
2004)  with  the  exception  that  radiative  fluxes  are  nonzero 
at  the  bottom  wherever  significant  radiation  can  penetrate 
to  that  depth.  In  this  situation,  the  radiation  reaching  the 
bottom  is  assumed  to  heat  the  bottom  layer  of  the  model 
and  also  provide  a  destabilizing  buoyancy  flux  that 
generates  turbulence  in  the  bottom  boundary  layer.  This 
bottom  buoyancy  flux  is  significant  only  in  very  shallow 


nearshore  regions.  Since  isopycnic  bottom  layers  are 
usually  much  thicker  that  the  bottom  boundary  layer 
(BBL)  in  the  deep  ocean,  use  of  the  BBL  parameterization 
is  optional  and  usually  invoked  in  coastal  ocean 
simulations  where  u  coordinates  provide  good  surface 
to  bottom  resolution  (Halliwell,  in  preparation).  When 
invoked,  BBL  mixing  is  only  implemented  at  grid  points 
where  at  least  one  vertical  coordinate  exists  within  the 
diagnosed  bottom  boundary  layer  thickness  above  the 
bottom.  Since  isopycnic  coordinates  migrate  to  resolve 
the  density  front  on  top  of  dense  overflows,  tests  are 
underway  to  determine  if  the  KPP  BBL  parameterization 
improves  the  representation  of  these  flows  (Section  2.3). 

The  three  vertical  diffusion  mixing  submodels  are 
capable  of  resolving  both  geostrophie  shear  and  ageos- 
trophic  wind-driven  shear  in  the  upper  ocean,  which  was 
not  possible  with  the  bulk  mixed  layer  of  MICOM. 
Halliwell  (2004)  demonstrated  this  by  forcing  HYCOM 
with  slowly  varying  monthly  climatological  forcing  in  a 
30-layer  Atlantic  Ocean  simulation  designed  so  that 
several  layers  were  available  to  resolve  the  surface  Ekman 
layer.  The  expected  Ekman  spiral  is  verified  at  two  model 
grid  points:  CRBN  in  the  Caribbean  Sea  and  NAC  in  the 
North  Atlantic  Current  (Fig.  3),  with  the  former  repre¬ 
senting  the  Trade  Wind  belt  and  the  latter  representing  the 
Westerly  Wind  belt.  Vector  velocity  relative  to  velocity  at 
the  Ekman  layer  base  (Fig.  3)  resemble  Ekman  spirals, 
indicating  that  although  geostrophic  velocity  shear  is 
present  in  all  of  the  velocity  profiles,  it  is  too  small  to  mask 
the  Ekman  spiral  structure.  Differences  among  the  spiral 
structures  are  associated  with  different  viscosity  profiles 
calculated  by  the  mixing  models  (not  shown).  The  central 
depths  of  the  reference  layers  can  be  used  as  a  proxy  for 
Ekman  layer  thickness  for  comparison  among  the  cases 
plotted  in  Fig.  3,  the  theoretical  Ekman  layer  thickness 
being  an  e-folding  scale.  At  point  CRBN,  the  MY  Ekman 
layer  is  thicker  than  the  KPP  and  GISS  Ekman  layers 
(67  m  versus  4 1  m  reference  layer  depth)  because  the  MY 
model  produces  larger  viscosity  coefficients.  The  NAC 
Ekman  layers  are  equally  thick  (78  m  reference  layer 
depth)  for  all  mixing  models  and  thicker  than  all  of  the 
CRBN  Ekman  layers.  Although  the  increasing  Coriolis 
parameter  acts  to  reduce  Ekman  layer  thickness  toward 
higher  latitudes,  the  larger  viscosity  present  at  point  NAC 
more  than  compensates  for  this  influence. 

To  illustrate  the  performance  of  two  vertical  mixing 
choices  in  continental  shelf  simulations,  zonal  sections  of 
temperature  and  vertical  viscosity  coefficient  across  the 
West  Florida  Shelf  at  27,55°N  are  presented  for  two 
HYCOM  experiments:  one  using  KPP  mixing  that  in¬ 
cludes  the  new  bottom  boundary  layer  parameterization 
and  one  using  MY  mixing  (Fig.  4).  Sections  are  shown  for 
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KPP  CRBN  VeL  Relative  to  Layer  9  KPP  HAC  VeL  Relative  to  Layer  1 3 


Fig.  3.  Winter  (mid -February)  velocity  vectors  (m/s)  at  two  model  grid  points:  CRBN  (left)  and  NAC  (right)  for  coarse  resolution  North  Atlantic 
HYCOM  simulations  (Halliwell,  2004)  using  KPP  (top),  GISS  (middle)  and  MY  (bottom)  mixing.  Vectors  are  shown  for  model  layers  located  above 
a  reference  layer  given  in  the  label  for  each  panel  and  chosen  by  inspection  to  reside  at  the  base  of  the  Ekman  layer  The  reference  layer  velocity  has 
been  subtracted  from  all  vectors  in  each  panel. 


March  22, 2002  when  there  was  a  strong  upwellmg  event 
in  the  presence  of  moderate  stratification*  The  magnitude 
and  distribution  of  vertical  viscosity  coefficients  is 
broadly  similar  between  the  KPP  and  MY  cases,  with 
distinct  surface  and  bottom  boundary  [ayers  present  over 
the  middle  and  outer  shelf.  Quantitative  differences  exist, 
with  weaker  turbulence  in  the  bottom  boundary  layer 
present  in  the  MY  case  compared  to  the  KPP  case. 


Stronger  nearbottom  stratification  in  the  MY  case  is 
associated  with  the  weaker  turbulence.  Offshore,  stronger 
turbulence  is  produced  in  the  surface  boundary  layer  by 
MY  compared  to  KPP.  This  conclusion  that  the  KPP  and 
MY  produce  qualitatively  similar,  but  quantitatively 
different  turbulence  patterns  has  also  been  reached  in 
idealized  shelf  simulations  using  POM  (Wijesekera  et  aL, 
2003)  and  ROMS  (Durski  et  aL,  2004). 
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Fig,  4.  Temperature  (left)  and  vertical  viscosity  coefficient  (right)  in  a  section  across  the  West  Florida  Shelf  from  two  simulations:  one  using  KPP 
mixing  with  the  new  bottom  boundary  layer  parameterization  (top)  and  one  using  MY  mixing  (bottom). 


2.3 .  Overflows 

A  proper  representation  of  overflow  waters  has  always 
been  challenging  for  OGCMs.  The  primary  reason  for  this 
difficulty  is  that  most  current  model  configurations  utilize 
horizontal  resolutions  that  cannot  explicitly  resolve  the 
complex  geometry  associated  with  most  overflows,  A 
further  challenge  in  the  modeling  of  overflows  is  that 
different  model  formulations  have  different  levels  of 
success  in  representing  them  (Gri  fries  et  aL,  2000).  Of 
particular  importance  is  the  model's  vertical  discretization 
(DYNAMO;  Willebrand  et  aL*  2001),  On  one  hand, 
terrain- folio  wing  (er)  coordinates  provide  the  ability  to 
concentrate  resolution  near  the  bottom  boundaries,  and 
hence  can  resolve  overflow  processes  quite  well  (Jung- 
ctaus  and  Mellon  2000),  provided  that  the  vertical  and 
horizontal  resolution  is  sufficiently  fine.  However,  the 
pressure-gradient  errors  associated  with  a  coordinates 
become  large  when  the  topography  is  steep.  Without  an 
explicit  representation  of  the  BBL,  z-  (or  p -)  coordinate 
models  tend  to  exhibit  unphysically  strong  entrainment  as 
gravity  currents  descend,  unless  both  the  vertical  reso¬ 
lution  is  fine  enough  to  resolve  the  BBL  thickness  (of 
order  tens  of  meters)  and  the  horizontal  resolution  can 
resolve  the  BBL  thickness  divided  by  the  slope  (of  order 
kilometers)  (Winton  etaL,  1998);  resolutions  much  higher 
than  currently  computationally  affordable,  Isopycnic 
coordinate  models  (or  hybrid  coordinate  models  that  are 
essentially  isopycnic  at  depth  such  as  HYCOM),  on  the 
other  hand,  have  vertical  resolution  that  naturally  migrates 


to  the  density  front  atop  a  gravity  current  and  do  not 
require  a  deviation  from  the  underlying  model  framework 
to  capture  the  structure  of  the  gravity  current  (Hallberg, 
2000).  In  isopycnic  coordinates,  there  is  no  numerically 
induced  diapycnal  mixing  and  it  is  necessary  to  explicitly 
parameterize  the  amount  of  mixing  occurring  during 
entrainment.  For  example.  Fig.  5  shows  the  colder  fresher 
water  forming  over  the  shelf  in  the  Nordic  Seas.  It  spills 
over  the  Denmark  Strait  entraining  more  saline  Irminger 
Sea  water  in  a  1/12°  North  Atlantic  HYCOM.  The  default 
parameterization  used  for  the  ocean  interior  in  that 
simulation  is  based  on  the  original  KPP  without  BBL 
parameterization.  Comparison  to  high-order  non-hydro- 
static  spectral  element  simulations  and  to  the  laboratory 
experiment  of  Turner  (1986)  however  shows  that  the 
interior  KPP  parameterization  (primarily  shear  instability 
mixing  tuned  for  the  ocean  interior)  underestimates  the 
amount  of  mixing  that  is  needed  for  entrainment,  but 
that  the  KPP  bottom  boundary  layer  parameters  can  be 
calibrated  to  agree  with  the  high-order  simulations 
(Chang  et  al.,  2005),  although  it  is  not  known  if  a  unique 
calibration  can  be  found.  This  is  illustrated  in  Fig.  6  which 
shows  the  Mediterranean  outflow  representation  in  a  1/ 
12°  regional  configuration  of  HYCOM  using  the  original 
KPP  parameterization  and  one  modified  to  increase 
mixing  in  high  isopycnal  slope  regions  (courtesy  of  X. 
Xu).  With  the  modified  parameterization,  the  Mediterra¬ 
nean  outflow  finds  a  neutral  depth  comparable  to  that 
observed  (which  is  too  deep  with  the  original  KPP).  More 
work,  however,  is  needed  to  fully  develop  a  physically- 
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Fig,  5.  Example  of  entrainment  in  the  Denmark  Straits  overflow  along  3l°W  in  a  cl imato logically -forced  l/12&  North  Atlantic  HYCOM:  colder 
fresher  water  forms  over  the  shelf  in  the  Nordic  Seas  and  spills  over  the  Denmark  Strait  and  entrains  more  saline  Irminger  Sea  water  (top  panel: 
temperature;  bottom  panel:  salinity). 


based  parameterization  that  can  be  used  in  all  regions  of 
high  vertical  shear. 

3.  The  prototype  Atlantic  Ocean  data  assimilative 
system 

While  HYCOM  is  a  highly  sophisticated  model, 
including  a  large  suite  of  physical  processes  and  incor¬ 
porating  numerical  techniques  that  are  optimal  for 
dynamically  different  regions  of  the  ocean,  data  assi¬ 
milation  is  still  essential  for  ocean  prediction  (a)  because 
many  ocean  phenomena  are  due  to  flow  instabilities  and 
thus  are  not  a  deterministic  response  to  atmospheric 
forcing,  (b)  because  of  errors  in  the  atmospheric  forcing, 
and  (c)  because  of  ocean  model  imperfections,  including 
limitations  in  resolution.  One  large  body  of  data  is 
obtained  remotely  from  instruments  aboard  satellites. 
They  provide  substantial  information  about  the  ocean's 
space-time  variability  at  the  surface,  but  they  are 
insufficient  by  themselves  for  specifying  the  subsurface 
variability.  Another  significant  body  of  data  is  in  the 
form  of  vertical  profiles  from  XBTs,  CTDs,  and  pro¬ 
filing  floats  (e.g.,  ARGO).  While  these  are  too  sparse  to 


characterize  the  horizontal  variability,  they  provide 
valuable  information  about  the  vertical  stratification. 
Even  together,  these  data  sets  are  insufficient  to  deter¬ 
mine  the  state  of  the  ocean  completely,  so  it  is  necessary 
to  exploit  prior  knowledge  in  the  form  of  statistics 
determined  from  past  observations  as  well  as  our  un¬ 
derstanding  of  ocean  dynamics.  By  combining  all  of 
these  observations  through  data  assimilation  into  an 
ocean  model  it  is  possible  to  produce  a  dynamically 
consistent  depiction  of  the  ocean.  It  is  important  that  the 
ocean  model  component  of  the  forecast  system  has  skill 
in  hindcasting  and  predicting  the  ocean  features  of 
interest.  Then  the  model  can  act  as  an  efficient  dy¬ 
namical  interpolator  of  the  observations.  The  1/16°  near 
global  Navy  Layered  Ocean  Model  (NLOM)  is  an 
example  of  how  an  ocean  model  can  be  a  successful 
dynamical  interpolator  of  surface  information  in  the 
assimilation  of  satellite  altimetry  observations  (Smed- 
stad  et  ah,  2003),  Shriver  et  al.  (2006-this  issue)  show 
that  the  1/32°  version  of  NLOM  is  an  even  better  dy¬ 
namical  interpolator. 

Performance  of  HYCOM  in  the  North  and  Equatorial 
Atlantic  has  been  documented  by  Chassignel  et  ah 
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Fig.  6.  Mediterranean  overflow  at  36°N  in  a  1/12°  regional  configuration  (which  includes  most  of  the  Gulf  of  Cadiz,  part  of  the  Hastcm  North 
Atlantic,  and  a  small  pan  of  the  Mediterranean  Sea);  left  panels,  original  KPP;  middle  panels  climatology;  right  panels,  modified  KPP.  Top  panels, 
temperature;  bottom  panels,  salinity  (courtesy  of  X.  Xu). 


(2003)  within  the  framework  of  the  Community 
Modeling  Experiment  (CME).  The  near  real-time  1/12° 
(—7  km  mid-latitude  resolution)  HYCOM  Atlantic 
Ocean  data  assimilative  system  (http://hycom.rs mas. 
miami.edu/oeean_prediction,html)  spans  from  28°S  to 
70°N,  including  the  Mediterranean  Sea  and  has  been 
running  since  July  2002.  The  vertical  resolution  consists 
of  26  hybrid  layers,  with  the  top  layer  typically  at  its 
minimum  thickness  of  3  m  (Le,,  in  fixed  coordinate  mode 
to  provide  near  surface  values).  In  coastal  waters,  there 
are  up  to  15  ^-levels  and  the  coastline  is  at  the  10  m 
isobath.  The  northern  and  southern  boundaries  are 


treated  as  closed,  but  are  outfitted  with  3°  buffer  zones 
in  which  temperature,  salinity,  and  pressure  are  linearly 
relaxed  toward  their  seasonally  varying  climatological 
values.  Three-hourly  wind  and  daily  thermal  forcing 
(interpolated  to  3  h)  are  presently  provided  by  the 
FNMOC  Navy  Operational  Global  Atmospheric  Predic¬ 
tion  System  (NOGAPS),  available  from  NAVOCEANO 
and  the  U.S,  GODAE  data  server  in  Monterey.  The 
HYCOM  data  assimilative  system  uses  surface  w'ind 
stress,  air  temperature,  and  specific  humidity  (from  dew 
point  temperature  and  sea  level  pressure)  in  addition  to 
shortwave  and  longwave  radiation.  Surface  heat  flux  is 
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Fig.  7.  (a)  The  sea  surface  height  from  the  1/1 2°  Atlantic  HYCOM  in  the  Gulf  Stream  region  on  April  1 3,  2005.  (b)  The  sea  surface  height  from  the 
1/12°  Atlantic  HYCOM  in  the  Gulf  of  Mexico  region  on  April  13,  2005.  Overlain  is  an  independent  frontal  analysis  of  high  resolution  MCSST 
observations  performed  at  NAVOCEANO.  The  frontal  position  is  marked  in  black  if  the  observations  are  more  than  4  days  old.  There  is  a  very  good 
agreement  between  the  model  frontal  location  and  the  one  determined  from  the  MCSST  frontal  analysis. 


calculated  via  a  bulk  parameterization  from  the 
NOGAPS  fields  and  model  SST.  The  model  was  spun- 
up  for  3  years  starting  from  the  end  of  a  previous  1/12° 
North  Atlantic  simulation  performed  with  the  M1COM 
(Chassignet  and  Garraffo,  2001). 

Mostly  because  of  its  simplicity,  robustness,  and  low 
computational  costs,  operational  ocean  prediction  sys¬ 
tems  around  the  world  (NLOM,  MERCATOR,  FOAM, 
etc.)  are  presently  using  Optimal  Interpolation  (OI) 
based  assimilation  techniques.  For  our  current  1/12° 
North  Atlantic  HYCOM  ocean  forecasting  system,  we 
have  adopted  a  similar  approach  by  selecting  an  OI 
technique  with  Cooper  and  Haines  ( 1 996)  for  downward 
projection  of  SSH  from  aJtimetiy.  The  basic  principle 
behind  Cooper  and  Haines  (1996)  is  conservation  of  the 
subsurface  potential  vorticity  during  the  assimilation 
time  step.  The  implementation  in  HYCOM  follows  that 
of  Hoang  et  a!.  (1997)  for  M1COM,  i,e,,  (1)  there  is  a 
reference  isopycnal  in  the  deep  ocean  below  which  the 
Montgomery  potential  is  unmodified  by  the  vertical 
projection  of  the  SSH  correction,  (2)  the  subsurface 
potential  vorticity  is  conserved  by  displacing  homo¬ 
genously  the  coordinate  interfaces  below  the  mixed 
layer  and  above  the  reference  isopycnal,  and  (3)  the 
geostrophic  relation  is  used  to  calculate  the  velocity 
increments  (except  in  the  equatorial  band  where  no 
velocity  increments  are  computed). 

Real-time  satellite  altimeter  data  (Geosat-Fol  tow-On 
(GFO),  ENVISAT,  and  Jason- 1 )  are  provided  via  the 


Altimeter  Data  Fusion  Center  (ADFC)  at  NAVOCEANO 
to  generate  the  two-dimensional  Modular  Ocean  Data 
Assimilation  System  (MODAS)  SSH  (1/4°)  analysis 
(Fox  et  al,,  2002)  that  is  assimilated  daily.  The  MODAS 
analysis  is  an  OI  technique  which  is  using  a  complex 
covariance  function  that  includes  spatially  varying 
length  and  time  scales  as  well  as  propagation  terms 
derived  from  many  years  of  altimetry  (Jacobs  et  al, 
2001).  Before  one  can  assimilate  the  SSH  anomalies 
determined  from  the  satellite  altimeter  data,  it  is  nec¬ 
essary  to  know  the  oceanic  mean  SSH  over  the  time 
period  of  the  altimeter  observations.  Unfortunately,  the 
geo  id  is  not  presently  known  accurately  on  the 
mesoscale.  Several  satellite  missions  are  either  under¬ 
way  or  planned  to  try  to  determine  a  more  accurate 
geoid,  but  until  the  measurements  become  accurate  to 
within  a  few  centimeters  on  scales  down  to  approxi¬ 
mately  30  km,  one  has  to  define  a  mean  oceanic  SSH.  At 
the  scales  of  interest  (tens  of  kilometers),  it  is  necessary 
to  have  the  mean  of  major  ocean  currents  and  associated 
SSH  fronts  sharply  defined.  This  is  not  feasible  from 
coarse  hydrographic  climatologies  (~  1  °  horizontal 
resolution)  and  the  approach  taken  by  the  HYCOM 
data  assimilative  system  has  been  to  use  a  model- 
derived  mean  SSH.  This  requires  a  fully  eddy-resolving 
ocean  model  which  is  consistent  with  hydrographic 
climatologies  and  with  fronts  in  the  correct  position. 
The  HYCOM-based  system  uses  the  model  mean  gen¬ 
erated  by  a  previous  1/12°  North  Atlantic  simulation 
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performed  with  MICOM  (see  Chassignet  and  GarrafTb, 
2001  for  a  discussion). 

The  model  sea  surface  temperature  is  relaxed  to  the 
daily  MOD  AS  1/8°  SSI  analysis  which  uses  daily 
Multi-Channel  Sea  Surface  Temperature  (MCSST)  data 
derived  from  the  5-channel  Advanced  Very  High  Reso¬ 
lution  Radiometers  (AVHRR)— globally  at  8,8  km 
resolution  and  at  2  km  in  selected  regions.  The  e-fold¬ 
ing  relaxation  time  is  a  function  of  the  mixed  layer  depth 
(30  days  x  mixed  layer  depth/20  m)  so  that  the  relaxation 
time  is  30  days  when  the  mixed  layer  depth  is  20  m  and 
300  days  when  the  mixed  layer  depth  is  200  m. 
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The  system  runs  once  a  week  every  Wednesday  and 
consists  of  a  1 0-day  hindcast  and  a  14-day  forecast.  The 
atmospheric  forcing  for  the  14-day  forecasts  gradually 
reverts  toward  climatology  after  5  days.  The  last  forecast 
record  is  weighted  with  the  contemporaneous  climato¬ 
logical  values  over  a  10  day  time  span.  Over  that  time,  a 
linearly  decreasing  (increasing)  weight  (1 -weight)  is 
used  for  the  forecast  (climatology).  During  the  forecast 
period,  the  SST  is  relaxed  toward  climatologically- 
corrected  persistence  of  the  nowcast  SST  with  a  relax¬ 
ation  time  scale  of  1/4  the  forecast  length  (i.e.,  1  day  for 
a  4-day  forecast).  The  impact  of  these  choices  is 
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Fig.  8.  The  sea  surface  height  (contour  interval  of  10  cm)  in  the  Gulf  of  Mexico  from  four  different  real-time  or  near  real-time  systems  overlain  on  the 
SeaWiFS  chlorophyll  concentrations  on  August  8,  2003.  Red/yellow  colors  indicate  areas  with  high  concentration,  while  the  darker  blue  color 
represents  areas  with  low  concentrations.  With  most  of  a  previously  detached  ring  reattached,  the  Loop  Current  is  elongated  and  extends  quite  far  to 
the  west  and  to  the  north.  Both  HYCOM  1/12°  and  the  global  NLOM  1/32°  do  a  good  job  at  capturing  the  full  northwestward  extent  of  the  Loop 
Current.  There  arc  some  small  differences  in  the  representation  of  the  recaptured  Loop  Current  nng;  the  ring  in  HYCOM  1/1 2°  is  still  showing  closed 
SSH  contours  with  a  signature  slightly  farther  north.  Both  HYCOM  1/12°  and  NLOM  1/32°  fail  to  correctly  place  the  eastward  frontal  position  of  the 
Loop  Current,  which  is  well  delineated  by  high  chlorophyll  in  the  observations.  In  the  global  NLOM  1/16°,  the  ring  did  not  remain  attached  and 
moved  westward  of  90*W.  The  Loop  Current  position  is  however  reasonably  well  represented  east  of  88°W.  In  the  1/24*  Navy  Coastal  Ocean  Model 
(NCOM)  configured  For  the  In  tra- Americas  Sea  (IAS),  the  Loop  Current  is  in  generally  the  right  location,  except  that  it  does  not  penetrate  far  enough 
northward  and  the  recaptured  ring  is  too  far  south  by  half  a  degree. 
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discussed  by  Smedstad  et  aL  (2003)  and  Shriver  et  al. 
(2006-this  issue), 

3.1.  Evaluation 

At  the  present  time,  evaluation  of  the  model  outputs 
relies  on  systematic  verification  of  key  parameters  and 
computation  of  statistical  indexes  by  reference  to  both 
climatological  and  real-time  data,  and,  in  a  delayed 
mode,  to  quality  controlled  observations.  The  accuracy 
of  data  assimilative  model  products  is  theoretically  a 
non-decreasing  function  of  the  amount  of  data  that  is 
assimilated.  A  degradation  caused  by  assimilation  gen¬ 
erally  indicates  inaccurate  assumptions  in  the  assimila¬ 
tion  scheme.  While  models  can  be  forced  to  agree  with 
observations  (e.g.,  by  replacing  equivalent  model  fields 
with  data),  improvements  with  respect  to  independent 
observations  are  not  trivial.  An  assessment  of  model 
improvement  (or  lack  of  degradation)  with  respect  to 
unassimilated,  independent  measurements  is  therefore 
an  effective  means  of  assessing  the  performance  of  an 
assimilation  system.  Variances  of  these  model-data  dif¬ 
ferences  serve  as  common  measures  of  the  estimation 
accuracy.  For  the  evaluation  of  flow  accuracy  and  water 
mass  characteristics,  we  follow  the  guidelines  put 
forward  by  the  international  GODAE  metrics  group 
(Le  Provost  et  al.,  2002)  as  well  as  the  validation  tests 
commonly  used  at  the  operational  centers  before  official 
transition  to  operational  use. 

Furthermore,  within  GODAE,  the  Atlantic  Ocean  has 
been  chosen  as  a  pilot  project  for  an  inter-comparison  of 


different  ocean  data  assimilative  systems,  including  tests 
and  evaluation  of  the  inter-comparison  process.  The 
Atlantic  was  chosen  because  of  the  developmental  status 
of  the  required  components  of  an  ocean  forecasting 
system:  already  well  instrumented,  large  number  of 
available  models,  high  user  interest.  The  comparison 
exercise  took  place  within  the  framework  of  MERSEA 
(Marine  EnviRonment  and  Security  for  the  European 
Area)  funded  by  the  European  Union  and  consisted  of 
comparing  similar  diagnostics  and  fields  from  corres¬ 
ponding  realizations  of  several  systems  (MERCATOR, 
FOAM,  MFS,  TOPAZ,  and  HYCOM)  (Crosnier  and  Le 
Provost,  2006-this  issue).  In  the  remainder  of  this  section, 
we  provide  evaluation  examples  for  the  HYCOM 
Atlantic  forecasting  system  that  differ  from  the  ones 
discussed  by  Crosnier  and  Le  Provost  (2006-this  issue). 

Examples  of  assessments  of  the  models'  ability  to 
represent  observed  flow  features  can  be  seen  in  Figs.  7a,  b 
and  8.  These  tests  qualitatively  evaluate  model  analyses 
against  alternate,  unassimilated  observations  of  flow 
features  in  regions  of  interest.  Fig.  7a  and  b  show  the 
model  SSH  hindcast  on  April  13,  2005  for  the  Gulf 
Stream  region  and  the  Gulf  of  Mexico  region.  Overlain 
on  the  SSH  is  an  independent  frontal  analysis  of  SST 
data.  Close  examination  of  the  figures  on  this  date  and 
others  (see  HYCOM  web  site  for  movies)  shows  overall  a 
very  good  agreement  between  the  model  frontal  struc¬ 
tures  and  the  independent  SST  observations.  Compar¬ 
isons  of  surface  height  and  temperature  with  ocean  color 
imagery  can  also  at  times  provide  clear  and  dramatic 
qualitative  model  assessment  (Chassignet  et  al„  2006). 


Sea  Surface  Temperaiure  at:  28  51  nQ84.5lw 


Sea  Surface  Temperature  at:  29nQ79w 
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Fig.  9.  Comparison  between  buoy  observations  of  SST  (blue),  the  1/12°  North  Atlantic  data  assimilative  system  (red),  and  MODAS  (black)  at 
(a)  28.5 1  °N,  84,5  i°W  and  (b)  29°N,  79CW. 
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Fig,  8  illustrates  the  value  of  SeaWiFS  ocean  color 
imagery  in  assessing  the  ability  of  four  ocean  prediction 
systems  to  map  mesoscale  variability  in  the  Gulf  of 
Mexico  and  in  helping  to  diagnose  specific  strengths  and 
weaknesses  of  the  systems.  Dark  areas  (chlorophyll  poor) 
are  found  in  the  interior  of  the  Loop  Current  and  a  semi¬ 
detached  eddy  while  light  areas  (chlorophyll  rich)  such 
as  the  Mississippi  River  plume  outline  the  eastern  edge  of 
the  eddy  and  the  northern  Loop  Current,  The  results 
clearly  illustrate  nowcast  accuracy  differences  in  the 
positioning  of  the  Loop  Current  and  in  representing  the 
eddy  by  the  four  prediction  systems  (see  Fig.  8  and 
Chassignet  et  al.,  2006  for  details). 

In  order  to  evaluate  whether  the  models  are  producing 
acceptable  nowcasts  and  forecasts  of  sea  surface 
temperature,  the  near  real-time  system  is  also  routinely 
compared  to  unassimilated  buoy  time  series  of  SST. 
Fig.  9  shows  two  examples;  the  first  (Fig,  9a)  is  a  buoy  in 
the  northern  part  of  the  Gulf  of  Mexico,  while  the  second 
one  (Tig.  9b  )  is  a  buoy  off  the  east  coast  of  Florida  in  the 


Atlantic  Ocean.  Time  evolutions  of  the  model  hmdeast 
SST  compare  well  to  the  independent  buoy  measure¬ 
ments  and  to  the  MODAS  SST  (the  latter  is  not  sur¬ 
prising  since  the  MODAS  SSTs  are  assimilated). 

Since  the  North  Atlantic  system  assimilates  only 
surface  quantities  (SSH,  SST),  quantitative  comparison 
of  model  temperature  and  salinity  to  unassimilated  pro¬ 
file  data  from  XBTs,  CTDs,  and  ARGO  floats,  and 
moored  buoys  allow  us  to  assess  system  performance  in 
the  ocean  interior,  including  the  skill  in  projecting 
surface  data  downward.  Two  good  examples  of  profiles 
are  shown  in  Fig,  10,  one  from  an  ARGO  profile  at 
5A°S  and  6,9°  W  (Fig.  10a)  and  one  from  a  PI  RATA 
buoy  at  10°N  and  10°W  (Fig,  10b).  Model  temperature 
sections  can  also  be  compared  to  XBT  measurements  as 
shown  in  Fig,  I !  for  the  Marine  Environmental  Data 
Service  (MEDS)  data  set.  A  quantitative  assessment  can 
then  be  performed  by  looking  at  the  RMS  difference 
between  the  model  and  data  profiles.  Fig.  12  shows  that, 
with  only  surface  data  assimilation,  the  version  of  the 


Saliniy  (psu) 


Fig.  10,  (a)  Temperature  and  salinity  profiles  from  ARGO  at  5.4°S  and  6.9°W  compared  to  the  1/12°  near  real-time  Atlantic  system  on  January  15, 
2004.  The  ARGO  temperature  profile  is  in  black,  the  salinity  profile  in  red.  The  corresponding  model  profiles  arc  in  blue  (temperature)  and  cyan 
(salinity),  (b)  Temperature  profiles  from  PI  RATA  at  I0°N  and  10°Won  April  13, 2004,  The  PlRATA  profile  is  in  red,  the  corresponding  model  profile 
in  cyan,  the  climatological  profile  in  black  and  the  MQDAS-3D  profile  in  blue. 
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Fig.  11.  (a)  Temperature  section  along  line  A  from  the  1/12"  near  real-time  Atlantic  system,  { b )  corresponding  section  from  the  MEDS  data, 
(c)  temperature  section  along  line  B  from  the  1/12°  near  real-time  Atlantic  system,  (d)  corresponding  section  from  the  MEDS  data. 


Atlantic  HYCOM  prediction  system  presented  here  has 
overall  a  greater  nowcast  RMS  error  than  climatology  or 
MODAS-3D.  MODAS-3D  (Fox  et  al„  2002)  uses  the 
statistics  of  the  historical  hydrographic  database  to 
downward  project  the  same  MODAS  SSH  anomaly  and 


SST  analyses  assimilated  by  HYCOM,  indicating 
superior  performance  for  a  data-based  method  of  down¬ 
ward  projection  than  the  Cooper  and  Haines  (1996) 
technique  used  in  HYCOM,  at  least  in  this  application. 
This  is  also  indicative  of  the  drift  in  T  and  S  that 
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Fig.  1 2.  (a)  Statistics  for  September  2003  through  August  2004  between  the  1/1 2*  HYCOM  system  and  available  Marine  Environmental  Data  Serv  ice 
(MEDS)  profile  observations.  The  RMS  difference  between  the  MEDS  data  and  different  climatologies,  MODAS  climatology  (CLIM),  MODAS 
synthetics  (MODAS),  Levitus  et  al.  ( 1 994),  and  the  Generalized  Digital  Environmental  Model  (GDEM3)  is  also  shown,  (b)  Statistics  for  September 
2003  through  August  2004  between  the  1/12°  HYCOM  system  and  available  PI  RATA  profile  observations.  The  RMS  between  the  PI  RATA  data  and 
the  MODAS  climatology  (CLIM)  and  MODAS  synthetics  (MODAS)  is  also  shown. 


occurred  during  the  spin-up  (Crosnier  and  Le  Provost, 
2006-this  issue).  Model  velocity  eross-seclions  can  be 
evaluated  through  qualitative  and  quantitative  compar¬ 
isons  of  biases  when  data  are  available.  Two  examples 
of  mean  velocity  comparisons  are  provided:  one  in 
Fig.  13  for  a  cross-equatorial  section  at  35°W  and  one  in 
Fig.  14  for  a  section  across  the  Yucatan  channel.  When 
observations  are  available,  transport  time  series  provide 
an  excellent  measure  of  the  model’s  ability  to  represent 
daily  to  seasonal  variability  (see  example  shown  in 
Fig.  1 5  for  the  Florida  Straits). 

3.2,  Model  outputs 

The  near  real-time  North  Atlantic  basin  model 
outputs  are  made  available  to  the  community  at  large 
within  24  h  via  the  Miami  Live  Access  Server  (LAS) 
(http://hycom.ramas.miami.edu/las).  Specifically,  the 
LAS  supports  model-data  and  model-model  compar¬ 
isons;  provides  HYCOM  subsets  to  coastal  or  regional 


nowcast/forecast  partners  as  boundary  conditions,  and 
increases  the  usability  of  HYCOM  results  by  “applica¬ 
tion  providers” 

4,  Boundary  conditions  for  regional  models 

An  important  attribute  of  the  data  assimilative  HYCOM 
simulations  is  the  capability  to  provide  boundary  condi¬ 
tions  to  regional  and  coastal  models.  The  chosen 
horizontal  and  vertical  resolution  for  the  forecasting 
system  marginally  resolves  the  coastal  ocean  (7  km  at 
mid-latitudes,  with  up  to  15  terrain-following  (a) 
coordinates  over  the  shelf),  but  is  an  excellent  starting 
point  for  even  higher  resolution  coastal  ocean  prediction 
systems.  To  increase  the  predictability  of  coastal  regimes, 
several  partners  within  the  HYCOM  consortium  are 
therefore  developing  and  evaluating  boundary  conditions 
for  coastal  prediction  models  based  on  the  HYCOM  data 
assimilative  system  outputs.  The  inner  nested  models  may 
or  may  not  be  HYCOM,  so  the  coupling  of  the  global  and 


76 


E.P.  Chassignet  et  aL  /  Journal  of  Marine  Systems  65  (2007)  60-83 


500 


1000 


1500 


2000 


2500 


3000 


1500 


4000 


4500 


5000 


EUC 

aesv 


SEC 

*2.0Sv 


NBC 

-6.GSv 

0 


SEC 
-7  0Sv 


Fig,  13.  (a)  Vertical  section  of  the  mean  velocity  across  the  Equator  at  35DW  from  5° S  to  5°N  from  the  1/1 2°  Atlantic  system  over  the  time  period  of 
September  2003  through  August  2004.  (b)  Observations  of  transports  across  the  Equator  at  35°W  (Schott  et  al.,  2003). 


coastal  models  needs  to  be  able  to  handle  unlike  vertical 
grids.  Coupling  HYCOM  to  HYCOM  is  now  routine  via 
one-way  nesting.  Outer  model  fields  are  interpolated  to 
the  horizontal  mesh  of  the  nested  model  throughout  the 


S7W  86W  85W 


entire  time  interval  of  the  nested  model  simulation  at  a 
time  interval  specified  by  the  user,  typically  once  per  day 
in  our  evaluations  to  date,  and  stored  in  HYCOM  archive 
format  Layers  can  be  added  to  these  archive  files  to 
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Fig.  14,  (a)  Vertical  section  of  the  mean  velocity  across  the  Yucatan  Channel  from  the  1/12°  Atlantic  system  over  the  time  period  of  September  2003 
through  August  2004.  (b)  Observ  ations  of  velocities  from  Abascal  et  al.  (2003), 
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Fig.  15.  The  transport  in  the  Florida  Current  at  27°N  from  the  1/12° 
Atlantic  near  real-time  system  is  in  black.  Observed  transport 
variations  in  the  Florida  Current  are  being  monitored  by  measuring 
the  cross -stream  voltages  using  an  undersea  cable  between  Florida  and 
the  Bahamas.  Daily  transport  data  arc  available  from  March  1982  to 
October  1998,  and  from  March  2000  onward  (Baringer  and  Larsen, 
2001).  Observations  from  the  cable  data  are  shown  in  red. 

increase  the  vertical  resolution  of  the  nested  model  and 
insure  that  there  is  sufficient  vertical  resolution  to  resolve 
the  bottom  boundary  layer.  The  nested  model  is  initialized 
from  the  first  archive  file  and  the  entire  set  of  archives 
provides  boundary  conditions  during  the  nested  run,  in¬ 
suring  consistency  between  initial  and  boundary  condi¬ 
tions.  This  procedure  has  proven  to  be  very  robust.  Nested 
Gulf  of  California  simulations  (Zamudio  et  al.,  submitted 
for  publication)  were  used  to  investigate  the  ability  of  the 


existing  HYCOM  nesting  capability  to  allow  accurate 
passage  of  a  coastally  trapped  wave  generated  by  Hurri¬ 
cane  Juliette  (Zamudio  et  aL,  2002).  In  Zamudio  et  al. 
(submitted  for  publication),  nesting  parameters  such  as 
updating  frequency,  e-folding  time,  and  buffer  zone  width 
were  varied  and  the  results  were  compared  to  coastal  tide 
gauge  stations.  Coupling  HYCOM  to  other  models,  such 
as  the  Navy  Coastal  Ocean  Model  (NCOM)  has  already 
been  demonstrated,  while  coupling  of  HYCOM  to 
unstructured  grid/finite  element  models  is  in  progress. 

In  the  remainder  of  this  section,  results  from  a  20-  layer 
1/25°  horizontal  resolution  Gulf  of  Mexico  HYCOM 
configuration  nested  within  a  non-assimilative  1999- 
2000  North  Atlantic  HYCOM  are  presented  and  dis¬ 
cussed.  The  model  domain  used  in  the  1/25°  nested  Gulf 
of  Mexico  HYCOM,  including  the  location  of  the  open 
boundary  conditions,  is  shown  in  Fig,  16.  Most  flow 
enters  the  domain  on  the  southern  and  southeastern  boun¬ 
daries  and  exits  through  the  Florida  Strait*  Currently  the 
nesting  of  HYCOM  to  HYCOM  is  one-way  (information 
is  only  passed  from  the  outer  grid  to  the  inner  grid)  and 
“off-line’;  meaning  that  the  nested  model  does  not  run 
concurrently  with  the  outer  model.  An  advantage  of  this 
approach  is  that  the  nested  region  does  not  need  to  be 
known  in  advance,  but  a  disadvantage  is  that  the  updating 
frequency  of  the  boundary  information  is  limited  by  how 
often  outer  grid  model  output  is  archived.  In  this  nested 
Gulf  of  Mexico  example,  the  method  of  characteristics 
(Browning  and  Kreiss,  1982,  1986)  is  used  for  the  baro- 
tropic  open  boundary  condition  on  velocity  and  pressure. 
At  the  open  boundaries,  20  grid  point-wide  “buffer^  (or 
boundary  relaxation)  zones  with  e-folding  times  of  0. 1  to 
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Fig.  16.  Bathymetry  used  in  the  1/25°  nested  Gulf  of  Mexico  simulation.  The  yellow  lines  indicate  the  locations  of  the  open  boundaries  that  are 
updated  from  a  1/12°  Atlantic  HYCOM  simulation. 
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10  days  (outer  to  inner  grid)  are  used  to  relax  the  baro- 
clinic  mode  temperature,  salinity,  pressure  and  velocity 
components  once  per  baroclinic  time  step  towards  a  non- 
assimilative  inter  annually  forced  1/12°  Atlantic  HYCOM 
solution  that  is  linearly  interpolated  in  time.  Although  the 
buffer  zone  is  located  on  the  fine  grid  mesh,  the  bottom 
topography  and  aforementioned  variables  are  constrained 
to  the  coarse  outer  grid  solution  and  thus  should  be  con¬ 
sidered  part  of  the  boundary  condition,  not  part  of  the  fine 
inner  grid  solution.  Concurrent  6-hourly  NOGAPS  was 
used  for  surface  forcing  in  both  the  nested  Gulf  of  Mexico 
model  and  the  interannually  forced  1/12°  Atlantic 
HYCOM  simulation. 

Compared  to  similar  1/12°  simulations,  the  1/25° 
simulation  shows  that  the  higher  resolution  results  in  more 
realistic  cyclonic  eddies  that  often  are  associated  with  the 
Loop  Current  and  Loop  Current  eddies  in  terms  of  eddy 


size,  speed,  population,  and  vertical  structure.  Fig,  17 
depicts  the  sea  surface  height  (top  panel )  and  sea  surface 
salinity  (bottom  panel)  on  June  1 3t  2000  (although  this  is 
a  non-assimilative  simulation  and  so  the  ocean  state  is  not 
representative  of  the  Gulf  of  Mexico  on  that  day).  At  this 
time,  the  Loop  Current  Extension  reaches  almost  28°N 
and  there  is  a  relatively  strong  cyclone  on  its  eastern  flank 
at  about  25°N.  This  cyclonic  eddy  plays  a  role  in  the  Loop 
Current  shedding  an  anticye  Ionic  ring  the  following 
month  (not  shown).  Several  of  the  cyclonic  eddies  travel 
along  the  Florida  Keys  and  then  exit  the  Gulf  of  Mexico 
through  the  Florida  Straits,  The  surface  salinity  shows  an 
area  of  relatively  fresh  water  just  north  of  the  Yucatan 
Peninsula  (88°W,  22°N).  This  is  an  area  of  prolific 
upwelling  during  June,  and  it  may  be  associated  with  the 
southward  Yucatan  Undercurrent  (Merino,  1997;  Ochoa 
et  al,  2001),  although  Cochrane  (1968,  1969)  suggested 
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Fig.  17.  Sea  surface  height  {top  panel)  and  salinity  (bottom  panel)  on  June  13,  2000  from  the  nested  1/25°  simulation. 
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that  bottom  friction  of  the  strong  Yucatan  Current  against 
the  slope  on  the  eastern  edge  of  the  Yucatan  Shelf  caused 
the  upwelling  instead. 

An  example  of  the  vertical  structure  of  the  anticyclonic 
Loop  Current  and  associated  cyclonic  eddies  on  June  1, 
2000  is  shown  in  Fig.  18.  At  this  time  the  Loop  Current 
extends  to  about  28°N  and  has  cyclones  on  both  the 
western  and  eastern  flanks.  A  zonal  cross-section  of  the 
south  to  north  (v)  component  of  velocity  and  the  salinity 
are  also  shown  in  Fig.  1 8.  The  eddy  on  the  eastern  flank  is 
propagating  along  the  continental  slope  but  in  this  case  the 
core  of  the  eddy  does  not  propagate  onto  the  shelf.  The 
salinity  depicts  a  subsurface  salinity  maximum  in  the 
center  of  the  anticyclone,  consistent  with  observations.  In 
all  of  the  nested  Gulf  of  Mexico  simulations,  realistic  flow 
through  the  southern  boundary  (about  4°  south  of  the 
Yucatan  Channel)  from  the  1/12°  Atlantic  HYCOM  is 
critical  for  realistic  Loop  Current  Eddy  shedding  in  the 
Gulf  of  Mexico.  In  particular,  the  flow  through  the 
Yucatan  Channel  needs  to  be  surface  intensified  on  the 
western  side  of  the  channel  and  have  a  mean  volume 
transport  of  about  28  Sverdrups  (Sv).  Although  recent 
measurements  (Sheinbaum  et  aL  2002)  suggest  the  value 
may  be  lower  than  this,  28  Sv  is  consistent  with  the  more 
extensively  observed  transport  through  the  Florida  Straits 
at  27°N  (Baringer  and  Larsen,  2001;  Johns  et  al,  2002). 

5,  Outlook 

The  long  term  goal  of  the  HYCOM  consortium  is 
an  eddy-resolving,  fully  global  ocean  prediction 


system  with  data  assimilation  to  be  transitioned  to 
the  ITS.  Naval  Oceanographic  Office  at  1/12° 
equatorial  (—7  km  mid-latitude)  resolution  in  2007 
and  1/25°  resolution  by  2011.  The  present  near  real¬ 
time  system  as  described  in  this  paper  is  a  first  step 
towards  the  fully  global  1/12°  HYCOM  data  assimi¬ 
lative  system.  Development  of  the  global  system  is 
presently  taking  place  and  includes  model  develop¬ 
ment,  data  assimilation,  and  ice  model  embedment. 
The  model  configuration  is  fully  global  with  the  Los 
Alamos  CICE  ice  model  embedded  and  will  run  at 
three  resolutions:  -60  km,  -20  km  and  —7  km  at 
mid-latitudes.  The  size  of  the  problem  makes  it  very 
difficult  to  use  sophisticated  assimilation  techniques 
since  some  of  these  methods  can  increase  the  cost  of 
running  the  model  by  a  factor  of  100.  The  strategy 
adopted  by  the  consortium  is  to  start  with  a  simple  data 
assimilation  approach  such  as  the  Cooper-Haincs 
technique  described  in  Section  4,  and  then  gradually 
increase  its  complexity.  Several  of  these  more 
sophisticated  data  assimilation  techniques  are  already 
in  place  and  are  in  the  process  of  being  evaluated. 
These  techniques  are,  ordered  by  degree  of  sophisti¬ 
cation,  the  NRL  Coupled  Ocean  Data  Assimilation 
(NCODA),  the  Singular  Evolutive  Extended  Kalman 
(SEEK)  filter,  the  Reduced  Order  Information  Filter 
(ROIF),  the  Reduced  Order  Adaptive  Filter  (ROAF) 
(including  adjoint),  the  Ensemble  Kalman  Filter 
(EnKF),  and  the  4D-VAR  Representer  method.  This 
does  not  mean  that  all  these  techniques  will  be  used 
operationally:  the  NCODA  and  SEEK  techniques  are 


Fig.  18.  Northeast  Guff  of  Mexico  zoom-in  of  sea  surface  height  (color),  surface  currents  (vectors)  and  bottom  topography  ("black  line  contours) 
(left  panel)  and  cross-sections  of  ^--component  of  velocity  (top  right  )  and  salinity  (bottom  light)  from  the  nested  1/25°  simulation  on  June  1 ,  20(H). 
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presently  being  considered  as  the  next  generation  data 
assimilation  to  be  used  in  the  near  real-time  system. 
The  remaining  techniques,  because  of  their  cost,  are 
being  evaluated  mostly  within  specific  limited  areas  of 
high  interest  or  coastal  HYCOM  configurations. 

The  NCODA  is  an  oceanographic  version  of  the 
multivariate  optimum  interpolation  (MVOI)  technique 
widely  used  in  operational  atmospheric  forecasting 
systems.  A  description  of  the  MVOI  technique  can  be 
found  in  Daley  (1991).  The  ocean  analysis  variables  in 
NCODA  are  temperature,  salinity,  geopotential  (dy¬ 
namic  height),  and  velocity.  The  horizontal  correlations 
are  multivariate  in  geopotential  and  velocity,  thereby 
permitting  adjustments  to  the  mass  field  to  be  correlated 
with  adjustments  to  the  flow  field.  NCODA  assimilates 
all  available  operational  sources  of  ocean  observations. 
This  includes  along  track  satellite  altimeter  observa¬ 
tions,  MCSST  and  in  situ  observations  of  SST  and  SSS, 
subsurface  temperature  and  salinity  profiles  from  BT’s 
and  profiling  floats,  and  sea  ice  concentration. 

Both  the  SEEK  filter  (Pham  et  al.,  1998)  and  ROIF 
(Chin  et  al.,  1 999)  are  sequential  in  nature,  implying  that 
only  past  observations  can  influence  the  current  estimate 
of  the  ocean  state  and  are  especially  well  suited  for  large 
dimensional  problems.  The  ROIF  assumes  a  tangent 
linear  approximation  to  the  system  dynamics,  while  the 
SEEK  filter  can  use  the  non-linear  model  to  propagate 
the  error  statistics  forward  in  time  (Ballabrera-Poy  et  al., 
200 1 ).  For  both  schemes,  the  analysis  step  is  multivariate 
in  nature,  i.e.,  all  model  state  variables  are  modified  in  a 
consistent  manner  after  the  analysis  step.  In  the  SEEK 
filter,  the  dominant  eigenvectors  describing  the  model 
variability  can  be  used  to  specify  the  initial  background 
error  covariance  matrix  in  decomposed  form.  This  leads 
to  fully  three-dimensional,  multivariate  dynamically 
consistent  corrections  (see  Parent  et  al.,  submitted  for 
publication  for  an  application  of  the  SEEK  filter  to  the 
North  Atlantic  configuration  of  Section  4).  The  ROI F 
method  factors  the  covariance  functions  into  horizontal 
and  vertical  components  and  represents  the  correction 
field  implicitly,  using  techniques  transplanted  from  sta¬ 
tistical  mechanics  (Gaussian  Markov  Random  Field). 
The  implicit  technique  tends  to  allow  a  highly  efficient 
way  to  represent  smaller  scale  dynamic  modes.  The 
reduced  order  aspect  of  ROIF  refers  to  the  fact  that  the 
information  matrix  is  approximated  as  a  banded  matrix. 
This  allows  more  realistic  tails  for  the  correlation  func¬ 
tions  than  similarly  approximating  the  error  covariance 
matrix. 

Finally,  another  Atlantic  configuration  is  under  de¬ 
velopment  to  form  the  backbone  of  the  NOAA/NCEP/ 
MMAB  North  Atlantic  Ocean  Forecast  System  (NAOFS). 


It  mostly  differs  from  the  system  described  in  Section  4  in 
two  ways:  (a)  different  horizontal  grid  and  (b)  NCEP- 
based  wind  and  thermal  forcing.  By  taking  advantage  of 
the  general  orthogonal  curvilinear  grid  in  HYCOM,  the 
NOAA/NCEP  group  is  using  a  configuration  which,  for 
the  same  number  of  grid  points  as  in  the  regular  Mercator 
projection  used  in  the  present  system,  has  finer  resolution 
in  the  western  and  northern  portions  of  the  basin  and 
on  shelves  (3-7  km)  in  order  to  provide  higher  reso¬ 
lution  along  the  US.  coast  than  toward  the  east  and 
southeast  (7-13  km). 
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Appendix  A.  Controlling  minimum  layer  thickness 
in  the  hybrid  coordinate  generator 


The  HYCOM  vertical  grid  is  controlled  by  reference 
isopycnals  and  the  minimum  thickness  permitted  for 
each  layer  k : 


dk  =  < 
where 

f  max[4*,min(4*,4)] 

l  4* 

(lsJfcSJV,)! 
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(3) 

4 =D/NS,  D  is  water  depth,  and  Ns  is  the  number  of  layers 
below  the  surface  permitted  to  transition  to  er  coordinates 
in  shallow  water.  To  estimate  4*  and  <$n,  the  minimum 
thicknesses  of  layer  1  (4i  and  4i),  the  largest  permitted 
minimum  thicknesses  (<5|tia*  and  ^“),  and  the  expansion/ 
contraction  factors  (ai  and  a2)  must  be  specified.  The 
expansion/contraction  factors  are  often  chosen  to  be 
greater  than  1  to  provide  the  highest  resolution  near  the 
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surface.  In  this  case,  the  minimum  thicknesses  will  increase 
with  depth  until  the  largest  permitted  values  are  reached, 
and  then  remain  constant  with  depth.  It  is  also  necessary 
to  identify  the  uppermost  model  layer  that  is  isopycnic 
(Ay  and  to  specify  an  interior  minimum  thickness  (Am,). 
Thickness  A  |*  governs  the  open  ocean  transition  between  P 
and  isopycnic  coordinates  by  maintaining  the  minimum 
thickness  of  the  nearsurface  P-coordinate  layers.  In  the 
isopycnic  coordinate  layers  beneath  the  P  domain,  mini¬ 
mum  thicknesses  are  specified  by  Am„  which  is  typically 
much  smaller  than  the  minimum  thicknesses  above  so  that 
sharp  pycnoctines  can  form  in  the  isopycnic  interior,  For 
the  present  experiments,  is  set  to  I  m. 

If  is  chosen  to  be  zero,  then  A*=A j*  everywhere  and 
the  coastal  transition  to  a  coordinates  is  not  implemented. 
If  Ns  is  nonzero,  then  the  coastal  transition  is  implemented 
in  model  layers  1  <k<Ns.  In  these  layers,  the  transition  to 
a  coordinates  begins  where  the  water  depth  becomes 
sufficiently  shallow  to  make  Where  the  water 

depth  is  shallow  enough  to  make  <5,  <<$2*,  a  transition  back 
to  level  coordinates  occurs  to  prevent  layers  from  be¬ 
coming  too  thin.  Parameter  is  updated  during  model 
runs.  The  thickness  of  layer  1  is  always  restored  to  a 
constant  thickness  A,  as  given  by  (1). 
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